Diffusion and percolation in anisotropic random barrier models 
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An anisotropic random barrier model is presented, in which the transition probabilities in diflerent 
directions have different probability density functions. At low temperatures, the anisotropic long- 
time diffusion coefficients, obtained using an effective medium approximation, follow an Arrhenius 
temperature dependence, with the same activation energy for each direction. Such activation energy 
is related to the anisotropic percolation properties of the lattice, and can be analysed in terms 
of the critical percolation path approximation. The anisotropic effective medium approximation 
is shown to predict the correct percolation threshold for an anisotropic two-dimensional square 
lattice. In addition, results are compared with numerical simulations using a fast kinetic Monte 
Carlo algorithm. 

PACS numbers: 05.40.-a, 05.60.Cd, 66.30.-h 



C/3 



X3 
O 

o 



> 
o 

O 

o 



o 



X 
S3 



I. INTRODUCTION 

Diffusion in disordered media is an active field of re- 
search, due to its relevance in a wide variety of natural 
and industrial processes P, |3j • One of the traditional 
models for disorder is the random barrier model (RBM), 
which consists of equally energy minima separated by en- 
ergy barriers, the height of which is randomly distributed 
according to a given probability density function (PDF). 
In this model, a particle moves from one minimum to 
another by performing thermally activated jumps. 

In these systems, diffusion properties can be studied ei- 
ther in time or frequency variables. Several studies have 
been conducted of diffusion properties in the isotropic 
RBM both under unbiased |i, |E S Sll Q and biased 
[inljl]| conditions. Following a frequency analysis (Refs. 
0,13 and references therein), the system may be charac- 
terized by a zero-frequency diffusion constant D[s — 0), 
and a characteristic frequency s* , which marks the onset 
of frequency-dependent diffusion. D{s = 0) and s* fol- 
low Arrhenius laws with the same activation energy Ec- 
Analogously, from a time variable standpoint, it takes a 
time t* ^ s*~^ for a particle to reach a long-time dif- 
fusion regime in the RBM, characterized by a diffusion 
constant D{t oo) = D{u) = 0) (Sj. The activation 
energy depends on the percolation properties of the 
lattice and the PDF of the energy barriers. This depen- 
dence is simply achieved by the critical percolation path 
approximation (CPPA), as shown for isotropic problems 
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In view of the diversity of systems in which diffusion 
takes place, the anisotropic generalization of diffusion 
problems has attracted considerable attention in the last 
years, both under unbiased El El El El IH H 113 
and biased [2^ l23l [2^ conditions. A few examples of 
anisotropic systems are porous reservoir rocks |^ JJ^ ^5] , 
layered semiconducting compounds |2^, and supercon- 
ductor cuprates ■ When dealing with anisotropic con- 
ditions, diffusion properties are independently studied in 
the different relevant directions of the system. It was 
recently shown that, for a two-dimensional anisotropic 



bond percolation model, different activation energies are 
found in each direction . 

In the present paper, a two-dimensional unbiased dif- 
fusion process is studied with each direction characterized 
by a different continuous PDF. The paper is organized as 
follows. In Sec. |n]the anisotropic RBM is introduced. 
In Sec. mil the model is studied within the anisotropic 
effective medium approximation (EM A). These results 
are compared against Monte Carlo simulations, whose 
numerical details are given in Sec. IIVI Section [3 is de- 
voted to a description of the CPPA ideas in anisotropic 
conditions, and in Sec. IVII the concluding remarks of the 
present paper are presented. 



II. ANISOTROPIC RANDOM BARRIER 
MODEL 



Diffusion processes will be studied on a two- 
dimensional square lattice with static disorder. Energy 
barriers are chosen from a given PDF p [E) aX, t — Q and 
are kept constant during the diffusion process. Possible 
jumps are only allowed between nearest neighbors. Once 
the energy barrier Eij between sites i and j is selected, 
the transition rates uJij from site i to site j are determined 



following an Arrhenius law 



Z 



(1) 



where ojq is the constant jump rate, z = 4 is the coordi- 
nation number, and (3 = l/ksT is the inverse tempera- 
ture, with ks being the Boltzmann constant. The energy 
Eij characterizes the bond joining sites i and j, therefore 
Eij = Eji, and the forward {i — > j) and backward {j ^ i) 
jumps have the same transition rate. 

In order to introduce the anisotropic character of the 
system, the Eij energies are selected from different PDFs, 
depending on the orientation of the bond joining sites i 
and j . Let 1 and 2 be the main directions of the square 
lattice, the key idea is to introduce pi{Ei) and p2{E2) 
instead of a single PDF p{E). The model is characterized 
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FIG. 1: Exponential (upper panel) and uniform (lower 
panel) probability density functions. Solid lines represent the 
anisotropic cases {a — 2), and dashed lines the isotropic cases. 



by anisotropy a = ei/e2 and global mean energy e = (ei + 
£2)72. The mean energies in each direction, ei and e2, are 
thus represented by ei = 2ae/{a + l) and €2 — 2e/{a + l). 
In the present work, a constant value of e = 0.5eo is 
adopted, where cq sets the unit of energy, and the effects 
of having a ^ I are studied. Two different anisotropic 
distributions will be considered: a) an exponential PDF 
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El £ [0, od) , 
i?2e[0,cx3), (2) 



and b) a uniform PDF 
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El e [(l-5i)ei,(l + (5i)ei]. 
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P2{E2)^tJ—, £^2 e [(1- 52)^2, (l+52)e2], (3) 
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where Si and 62 serve to control different distribution 
widths in each direction. This uniform PDF represents 
the most general anisotropic extension of the isotropic 
uniform PDF used in Ref. [g to study diffusion in RBM. 
In the following, and for the sake of simplicity, the widths 
of the uniform PDF will be 61 — S2 — 0.5. Figure^shows 
the exponential and uniform PDFs for a = I and a ~ 2. 



III. ANISOTROPIC EMA: LOW 
TEMPERATURE PREDICTIONS 

The EMA self-consistent conditions provide a method 
for obtaining the diffusion coefficients for a given disor- 



dered medium. Usually, these equations must be numer- 
ically solved, except for some simple cases. It is showed 
in this section that for the low temperature limit, some 
analytical predictions may be obtained within the RBM. 



A. Self-consistent conditions 

Many authors [H IH US iH have derived the EMA 
that allows to obtain effective diffusion coefficients. The 
approach considers one impure bond of the disordered 
lattice as embedded in an effective medium, mimicking 
the average surroundings. By imposing the averaged 
fluctuations to be zero, the self-consistent condition is 
derived for the transition rate of the effective medium. 
For an hypercubic d-dimcnsional isotropic lattice in the 
long-time limit, this condition reads j^]: 



+ {d~l)a 



0, 



(4) 



where to is the transition rate of the impure bond dis- 
tributed according to the PDF i/(w), a is the transition 
rate of the effective medium, and the brackets denote 
average over the PDF Solving the self-consistent 

condition for a, the diffusion coefhcient is obtained as 
D — a a? , where a is the lattice constant. 

The anisotropic extension of such formalism, where 
there exist n different directions, leads to n coupled equa- 
tions that self-consistently solve for the n different dif- 
fusion coefficients. In a two-dimensional square lattice, 
and for the long-time limit, these conditions are ilSul*^ 







0, 



with 



Jmn — — arctan ^ / — 
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(5) 



(6) 



and m, n = 1,2 denoting the principal axes of the lat- 
tice. The effective transition rates Ci are related to the 
diffusion constants by Di = aiO^ . 

At high temperatures, the particle can easily overcome 
energy barriers, and eventually all diffusion constants 
approach the same value. In Fig. |21 the normalized 
diffusion coefficients at high temperatures for the two- 
dimensional square lattice with an exponential PDF are 
plotted as functions of temperature, both under isotropic 
and anisotropic conditions. Lines represent the solutions 
of the EMA self-consistent conditions Eqs. I@J and jSJ, 
and symbols correspond to numerical simulations (see 
Sec. IIVII . The figure shows that, for high temperatures, 
Di/a^LdQ 1/z. Analogous results are obtained using 
the uniform PDF. In the next Subsections the predictions 
of EMA for diffusion at low temperatures are considered. 
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FIG. 2: Diffusion coefficients at higii temperatures. Lines 
correspond to the EMA solution and symbols to numerical 
simulations. The isotropic case is represented by a dashed 
line and the anisotropic case a = 2 by continuous lines. 



In the EMA, the bond percolation threshold of the hy- 
percubic lattice is given by pf^^ = .^^^IHIlJ. 
Therefore, Eq. is a condition over Ec for each partic- 
ular PDF p{E), in terms of the percolation properties of 
the lattice. Indeed, combining this value of Ec with Eq. 
((HI, the EMA diffusion coefficient for isotropic d dimen- 
sional hypercubic lattices at low temperatures is given 

by 

It is worth noting that ~ dr^ is only exact for 

d — 2 32], therefore Eq. H12|l does not give the correct 
exponential behavior for d — 2,, and other approxima- 
tions, such as CPPA, should be considered |2J- 

C. Anisotropic two-dimensional case 



B. Isotropic case 

For an isotropic hypercubic lattice in d dimensions, Eq. 
Q may be written as 



{uj + {d- l)a - da) 
UJ+{d- l)cr 



= 0. 



(7) 



By introducing the explicit dependence on energy u) = 
Lo{E) given in Eq. JQ), and transforming the transition 
rate average over v{uj) to an energy average over p{E), 
Eq. (UJ becomes 
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u){E) + (d- l)cr 



P(-E) 



1 

da 



(8) 



For /3 — > oo the transition rate u;(E) varies extremely 
rapidly, due to its exponential dependence. Therefore it 
is possible to define an energy value Ec such that two 
possibilities arise: uj{E) ^ {d — l)a for E < Ec, or 
uj{E) ^ (d — l)a for E > Ec- The characteristic value 
Er can be therefore defined as 



u;{Ec) = (rf - l)a. 



(9) 



For values of < Ec, the left hand side of Eq. ^ 
vanishes. Alternatively, for E > Ec the value uj{E) in 
the left hand side of Eq. (jHJ may be ignored. Taking 
these conditions into account, and averaging over p{E), 
Eq. ^ for /3 — > oo becomes 



1 

da 



P{E) 
{d-l)a 



dE. 



(10) 



Or, equivalently, 



jp{E)dE^\. 



(11) 



For the anisotropic two-dimensional case, Eq. (0 
turns into two self-consistent conditions, with transition 
rate PDFs and 1^2 (W2), for each direction of the lat- 

tice. By introducing the energy dependence uji = uj{Ei) 
and 0J2 = w(i?2), the two self-consistent conditions read 



(^^i) + if 12 - 1) ^1 



,^(^2) + (/2l'-l) 'T2, 



PiiEi) 



P2{E2) 
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111 
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(13) 



Again, the PDFs change abruptly for f3 00, and 
a parameter Ec can be defined as in the isotropic case. 
However, Ec is expected to be characteristic of the un- 
derlying energy landscape, so the diffusion coefhcients in 
each direction are expected to be governed by a single 
Ec- For the anisotropic case, an energy Ec will be de- 
fined separating two limiting conditions simultaneously: 
u> (El) < (/f2^ - 1) ai and to {E2) < {f2i - l) 0-2 for 
El and E2 larger than Ec, and to (Ei) :$> (/j^^ — l) ai 
and uj {E2) ^ (/2i^ ~ 1) o'2 for Ei and E2 smaller than 
Ec- Thus, Ec must verify two simultaneous conditions: 



ojiEc) = {fi2'^i)<yi, 

^{Ec) - if 21 - l)f^2. 



(14) 



In this way, a set of equations analogous to Eq. pi(l are 
obtained, 



J pi[Ei)dEi ^ h2, 


j p2{E2)dE2 = /21. 



(15) 
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By adding Eqs. H15|) . using Eq. © and trigonometric 
relations, an expression is arrived at, 



Be 



Pi{Ei)dEi 



P2{E2)dE2 = 1, 



(16) 



which gives the activation energy Ec as a function of 
the anisotropy a. Moreover, by replacing the expressions 
in Eqs. H15|) for /12 and /21 in Eqs. H14|) . and solving 
for (Ti and (72, the corresponding anisotropic diffusion 
coefficients are obtained, 



/ pi{Ei)dEi 




Do 



Ec 

/ p2iE2)dE2 


Ec 

/ p2{E2)dE2 2 

_0 ^QQ ^^-/3Ec 

Ec r 



(17) 



/ pi{Ei)dEi 




The isotropic result Eq. H12I) is obviously recovered by 
setting pi = p2 • 

FiguresOland^show Arrhenius plots of the anisotropic 
diffusion coefficients at low temperatures, corresponding 
to the exponential and uniform PDFs, respectively. The 
solution of EM A self-consistent conditions Eqs. Q and 
(jnj are represented with dashed an continuous lines, for 
the isotropic and anisotropic a — 2 cases, respectively. 
Dotted lines represent the prediction of EMA for the low 
temperature limit, Eqs. (|17|l . Symbols are the results of 
numerical simulations, as described in the next section. 
These figures show that the diffusion coefficients in each 
direction follow Arrhenius laws with the same activation 
energy. Even though simulation data at lower tempera- 
tures are needed, the agreement with the low tempera- 
ture anisotropic diffusion coefficients is better for a uni- 
form PDF than for an exponential PDF. 



IV. NUMERICAL SIMULATIONS 
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FIG. 3: Arrhenius plot of the diffusion coefficients for an 
exponential PDF. The solution of the self-consistent EMA 
conditions is represented with a dashed line for the isotropic 
case, and with continuous lines for the anisotropic a = 2 cases. 
Dotted lines represent the analytical EMA predictions, Eqs. 
117II . for low temperatures. Symbols correspond to SMC and 
FKMC simulations, as indicated. 
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FIG. 4: Arrhenius plot of the diffusion coefficients for a uni- 
form PDF. The solution of the self-consistent EMA condi- 
tions is represented with a dashed line for the isotropic case, 
and with continuous lines for the anisotropic a — 2 cases. 
Dotted lines represent the analytical EMA predictions, Eqs. 
117II . for low temperatures. Symbols correspond to FKMC 
simulations. 



Monte Carlo simulations were performed to obtain 
long-time diffusion coefficients for comparison with 
anisotropic EMA predictions. The energy landscape was 
selected from the corresponding PDF at t = and kept 
fixed during the diffusion process. At t = a particle was 
assigned to a random initial site i. Different Monte Carlo 
algorithms may be used at this point and two possibilities 
were considered: standard Monte Carlo (SMC) ^ and a 
fast kinetic Monte Carlo (FKMC) 33] scheme. A brief 
description of these methods is given in the following. 

In SMC, the particle selects at random one of its near- 
est neighbors j and tries to overcome the barrier between 
them in a time unit. A random number ^ S (0, 1) is gen- 
erated such that if ^ < cuij the jump is effective, otherwise 



the particle stays at the initial site. In this process, one 
unit of time is used for every jump trial. Although SMC 
simulations proved to be very useful for studying diffu- 
sion processes, it was shown that it is not too appropriate 
for studying low temperature regimes 5, 8, 33.j. At low 
temperatures, the transition rates decrease exponentially 
with increasing /3, and the random number ^ is, mostly, 
orders of magnitude greater than the transition rates, 
making the number of effective jumps (displacements) 
very small and the long-time diffusion regime difficult to 
reach. 

In the FKMC [23| scheme, consider the particle in a 
site i on a lattice with its z nearest neighbors j (j = 
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1, z). The transition rates from i to j are denoted u>ij. 
The total transition rate uj.: from site i is defined as: 



(18) 



Instead of selecting the neighbor at random, as in SMC, 
a neighbor k is selected for an effective jump given that 



1 ''^^ 1 

— Y] < Ci < — 



j=k 



(19) 



where is randomly uniformly distributed in (0, 1). The 
time variable t is then increased in t' , where t' is chosen 
from an exponential distribution with mean waiting time 
Therefore, 



UJi 



(20) 



with ^2 randomly uniformly distributed in (0,1). This 
procedure is repeated from site k and so on. In the 
FKMC algorithm, each jump trial is effective, meaning 
that the particle always jumps to one of its neighbors, 
and the time elapsed in one jump is accordingly adjusted. 
Furthermore, the FKMC algorithm depends on the ratios 
ujii/ uJi and consequently it is not as /3 dependent as cjji 
[3^ . This algorithm allows to reach larger values of (3 in 
simulations of the diffusion process. 

Simulations were performed on 300 x 300 and 500 x 500 
sites square lattices for the SMC and FKMC, respec- 
tively, with periodic boundary conditions. For each algo- 
rithm, the mean square displacements on each direction 
(r^(i)) and (^K^)) were computed, averaging over be- 
tween 2000 and 5000 realizations of the random walks. 
The long-time diffusion coefficients were defined through 

{''1(2) (^)) ^ 2Z)i(2)t, and were obtained from the best 
linear fits to the long-time mean square displacements. 

In Figs. 01 and 01 numerical simulations and EM A re- 
sults are presented together. In Fig. O SMC simulations 
are plotted up to Pcq — 10 and some FKMC simulation 
points are shown for comparison. Both algorithms co- 
incide within the numerical precision. In Fig. ^ only 
FKMC results are presented up to a value /3eo — 30. 
Monte Carlo simulations do not completely reach the 
asymptotic low temperatures behavior. However, numer- 
ical simulations and EMA seem to agree very well in the 
accessible temperature range. 



At low temperatures the characteristic Arrhenius dif- 
fusion energy Ec can be related to the bond percolation 
threshold of the lattice. Consider a random walk on a 
realization of the disorder energy landscape at a very 
low temperature. In order to overcome a barrier with 
an energy E' , the particle spends a mean waiting time 
t' ~ exp(/3i?'). For short times, therefore, the particle 
can only move to sites which are connected by low energy 
barriers and is surrounded by a perimeter of higher en- 
ergy. Roughly, at time t' the particle might jump barriers 
with E < E' , and the probability to overcome this bar- 
riers is p{E)dE. For longer times, the particle could 
overcome the lowest barrier of the perimeter, and access a 
new region with a higher energy perimeter. These regions 
are non-compact in the sense that they may have inside 
barriers that belong to the perimeter barriers. Eventu- 
ally, there exists a particular barrier of height Ec , beyond 
which the particle gains accesses to the whole system, 
through the corresponding percolation path of energies 
E < Ec- Thus, for an isotropic medium, Ec is given by 



/ 



p{E)dE = pc, 



(21) 



where pc is the bond percolation threshold of the system 
(pc = 0.5 for the two dimensional isotropic square lattice 
|32 | ) . It has been shown that Ec is the highest energy bar- 
rier which the particle must overcome in order to gain full 
access to the percolation network. The long-time diffu- 
sion coefficient must therefore he D ^ exp{—f3Ec), which 
is indeed the observed behavior of isotropic diffusion at 
low temperatures 0, . 

The percolation threshold of a particular lattice, which 
is given by a point pc for isotropic percolation, becomes 
for anisotropic percolation a critical surface fHpi}) = 
|33 |. where {pi} denotes the set of relevant occupation 
probabilities. For example, the percolation function is: 
ip{p) = p — Pc for isotropic percolation, (p{pi,p2) = 
Pi + P2 — ^ for the square lattice, and ^{pi,P2,P3)_ = 
Pi + P2 + P3 — P1P2P3 — 1 for the triangular lattice [3J| . 
Furthermore, the critical surface implies a change in the 
morphology of the incipient percolation network. 

In the anisotropic RBM context, the occupation prob- 
abilities of a bond with energy barrier E, i.e. accessi- 
bility condition of the bond, is given by the probability 
of E being lower than the maximum accessible barrier. 
Therefore, the generalization of Eq. H21|l to anisotropic 
conditions becomes 



V. CRITICAL PERCOLATION PATH 
APPROXIMATION 

The idea of a percolation path governing diffusion at 
low temperatures was first developed in Ref. |0| and 
rigorously proved later In this section, this idea 

will be briefly summarized and extended to anisotropic 
conditions. 




= 0. 



(22) 



Note that there exists just one energy Ec, which is the 
same for all directions, and gives full access to the whole 
anisotropic percolation network. For the anisotropic 
RBM on a square lattice, Eq. H22|) becomes equal to the 
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FIG. 5: Dependence of the activation energy Ec on the 
anisotropic parameter a for the exponential and uniform 
PDFs studied here. 

EMA prediction, Eq. (HBJ. This means that EMA pre- 
dicts the correct critical percolation surface ip(jji,p2) = 
Pi +P2 — 1 for anisotropic bond percolation in the square 
lattice ^,'21]. 

Figure |31 shows the effect of anisotropy a on Ec{a) for 
the energy distributions studied in the present model, 
and predicted both by CPPA Eq. and EMA Eq. 

(|16|l . For the exponential PDF, the condition for Ec reads 
exp(— i?c/ei) + exp(— £'c/e2) = 1, which was numerically 
solved. For the uniform PDF Eq. gives Ec/eo — 

2a/(a + l)2. 

VI. CONCLUDING REMARKS 

In this paper, diffusion properties were studied using 
an anisotropic RBM, with emphasis on the low temper- 
ature behavior and on percolation properties. Two kind 
of PDFs were used to characterize different directions of 
the lattice, namely, exponential and uniform PDFs. The 
anisotropic EMA was used to calculate the long-time dif- 
fusion properties for all temperatures, derived from the 
numerical solutions of the self-consistent conditions ex- 
pressed in Eqs. Q. Furthermore, analytical expressions 



for the diffusion coefficients at low temperatures were ob- 
tained, Eq. (|17|l . which show that diffusion in different 
directions follows Arrhenius laws with a same activation 
energy Ec- This should be compared with the thermally 
activated diffusion in anisotropic bond percolation lat- 
tices, in which different activation energies are found for 
each direction 18]. In the present model, only one activa- 
tion energy is found due to the existence of an anisotropic 
percolation path of low energy barriers, which governs 
the diffusion process. Besides of giving the activation 
energy for diffusion, EMA predicts the exponential pref- 
actor for diffusion and it dependence with the anisotropy 
of the disordered system. 

The two Monte Carlo algorithms used here, namely, 
SMC and FKMC, show a very good agreement with the 
EMA predictions for the diffusion coefficients in the ac- 
cessible temperature range. For a more extensive com- 
parison with EMA, other algorithms should be used. 

In the present paper, a connection was established be- 
tween EMA and CPPA ideas, and EMA was shown to 
predict the correct activation energy for anisotropic dif- 
fusion in the square lattice. For other geometries and 
dimensions, it is expected that EMA will still predict an 
Arrhenius behavior, but with an activation energy that 
differs from that predicted by CPPA. This difference is 
due to the fact that CPPA uses the percolation threshold 
as a parameter, while EMA predicts its own percolation 
threshold. However, EMA is known to predict the correct 
percolation threshold o nly for the square lattice, even in 
anisotropic conditions pol l2l| . Concerning CPPA, cor- 
rections of the form become relevant for dimensions 
grater than two, but it is not clear which of both approx- 
imations, EMA or CPPA, give better results J} and a 
systematic comparison turns necessary. Additional work 
on this direction is now under progress. 
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